///usage:  ini ROOT528
///usage:  root -l -q  readtribn-data.root 'myFit.C(0,0,0)' ## Only bkg models+data
///usage:  root -l -q  readtribn-data.root 'myFit.C(0,0,0,"merge_list")' ## Only bkg models+data , merge some bkg templates
///usage:
///usage:  myFit.C knows 3 arguments: 0,0,0 -- 1) 0,1  corresponds to Btag OP THCP6/THCPT, 2) 0,1 un-/corrected templates 
///usage: 3) scenario 0,1 TMinuit/TFitter MINOS-like algo 


///usage: script can write all output to log file
///usage: you can extract intresing points by doing
///usage: sed -ne '/\+MYOUTPUT/,/\-MYOUTPUT/p' myFit.log



#include "TObjString.h"
#include "TList.h"
#include "TMath.h"
#include "TLegend.h"
#include "TGraph.h"
#include "TArrayD.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TROOT.h"
#include <iostream>
#include "TMinuit.h"
#include "TFitter.h"

#include "TMatrixDSym.h"
#include "TGraphAnalysis.h"

#include "TGraph.h"

#include <unistd.h>
#include <stdio.h>

//#define doLogFile ///usage
#define LOGFILE "myFit.log" ///usage

///usage ON/OFF
///usage you can comment ==>
///usage  sed -e 's///#define doHESSE/\/\/\//#define doHESSE/g' myFit.C > myFit.C.tmp; mv myFit.C.tmp myFit.C;
///usage or uncomment ==>
///usage  sed -e 's/\/\/\//#define doHESSE///#define doHESSE/g' myFit.C > myFit.C.tmp; mv myFit.C.tmp myFit.C; 
#define doMIGRAD ///usage
#define doHESSE ///usage
#define doMINOS ///usage

///usage:
///usage: to make the table of relevant info
///usage: sed -ne '/\+MYOUTPUT/,/\-MYOUTPUT/p' myFit.log | grep -E "=|VALS|MIGRAD|MINOS|RATIOS" | tr  '=' ';' | tr -d ' ' | sed -e 's/_TCHPTTrig0//g' -e 's/yield//g' -e 's/ratio//g'
///usage:

#define NUMSCAN 100 ///usage: number of Chi2 scanning point
#define NARROWSCAN ///usage: make scan narrower
//#define NARROWRANGE 0.5 ///usage: make scan narrower

///To merge histograms
#include "templateMerger.h"

int colors [] =
{
        kRed,
        kGreen,
        kMagenta,
        kYellow,
        kRed-2,
        kGreen-2,
        kMagenta-2,
        kYellow-2,
        kRed+3,
        kGreen+3,
        kMagenta+3,
        kYellow+3,


};

///Templates && data histogram properties
TH1 * dataHist=0;
double minVal=-1e30;
double maxVal=+1e30;
double dataYield=1e0;
int bins = 0;

///Flag to understand the status of the current code execution
bool _isOk=false;



////Containers
TList * arrayOfHistos=0; /// First element is always data histogram
TList * arrayOfNames =0;
//TList * arrayOfYields=0;
TList * arrayOfLegends=0;


/// we will use two additional index arrays (with the same size) carring info about
/// 1) flavor 0,1,2 (for "Q", "C", "B");  size of values == nfc
/// 2) Cat 0,1,2 ( for Bbb,bBb,bbB); size of values == ncateg
/// sizes of flavors and positions arrays must be the same

/// Here


///     Bbb bBb bbb
///     Cbb bCb bbC
///     Qbb bQb bbQ

/// is coded

/// Carries 'Flavor' index
 const int flavors [] = {
        2, //B
 //    2, //B
//      2, //B
//        1, //C
     1, //C
//      1,  //C
//        0, //Q
      0, //Q
//      0  //Q
        };

/// Carries 'Cat' index
 const int positions [] = {
	0, //Bbb
//      1, //bBb
//       2, //bbB
//	0, //Cbb
        1, //bCb
//      2, //bbC
//	0, //Qbb
       1, //bQb
    //  2  //bbQ
        };





/**
                Performs merging 
                _mother == mother template ( ie bbB)
                _daughters == TList * of daughter templates (TObjString* names)
                _mode == mode of operation: Merge && Normalize -- 1, Normalize && Merge && Normalize 
        It changes :
                * arrayOfHistos
                * arrayOfNames
                deletes  corresponded _daughters templates references


**/

bool templateMerge(TString _mother, TList *_daughters, int mode=1)
{

        if ( _mother.Length()==0 && _daughters == 0 && _daughters->IsEmpty()) {std::cout<<"templateMerge --> something wrong"<<std::endl;return false;}

TList *_remove_arrayOfHistos = new TList();
TList *_remove_arrayOfNames = new TList();


///Find mother
TH1 *  _motherTh1=0;
for (int i=0; i<arrayOfNames->GetSize();i++)
{
TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
if (_name.Contains(_mother)){ _motherTh1 =  (TH1*)arrayOfHistos->At(i); break;}
}


if (_motherTh1==0) {std::cout<<"templateMerge --> something wrong mother"<<std::endl;return false;}

///Find daughter templates
for (int i=0; i<arrayOfNames->GetSize();i++)
{

TH1 * tempHis = (TH1*)arrayOfHistos->At(i);

TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();
 for (int jj=0; jj<_daughters->GetSize();jj++)
{
        TString _nameDaughter =  ( (TObjString *) _daughters->At(jj))->GetString();
        if (_name.Contains(_nameDaughter)) {

	cout<<"I'm merging " << _name<<"\n";

        _remove_arrayOfNames->Add(arrayOfNames->At(i));
        _remove_arrayOfHistos->Add(arrayOfHistos->At(i));

        } ///if
} ///for


} ///for


///Merging daughters to mother && Deleting daughter from arrayOfNames arrayOfHistos
if (!templateMerger(_motherTh1,_remove_arrayOfHistos,mode)) {std::cout<<"templateMerge --> something wrong with merging"<<std::endl;return false;}


for (int i=0; i<_remove_arrayOfNames->GetSize();i++)
{
        arrayOfHistos->Remove( _remove_arrayOfHistos->At(i));
        delete  _remove_arrayOfHistos->At(i);
        arrayOfNames->Remove( _remove_arrayOfNames->At(i));
        delete _remove_arrayOfNames->At(i);
}


delete _remove_arrayOfHistos;
delete _remove_arrayOfNames;
return true;

}





void chi_square(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  //calculate chisquare
  double chisq = 0;

TH1 * _DataHist = (TH1*)  arrayOfHistos->At(0);
TH1F * _combinedHist = new TH1F("combined_hist","combined_hist",_DataHist->GetNbinsX(),_DataHist->GetXaxis()->GetXmin(),_DataHist->GetXaxis()->GetXmax());


///Combine all templates in one 
for (int k=1;k<arrayOfHistos->GetSize();k++) {

//std::cout<<"A'm adding "<<arrayOfNames->At(k)->GetName()<<"with par "<<par[k-1]<<"\n";

_combinedHist->Add((TH1 * ) arrayOfHistos->At(k),par[k-1]);



}

for (int k=1;k<=_DataHist->GetNbinsX();k++) {

double _error =  _DataHist->GetBinError(k);

double delta=0;

if (_error>0) delta  = ( _DataHist->GetBinContent(k)-_combinedHist->GetBinContent(k))/_error;
chisq += delta*delta;

}
delete _combinedHist;

f=chisq;

  return;
}


/**
        Parameters of the function:
        _novl = 0,1 (TCHPT, TCHP6), Btag operation point of mass templates
        _corr= 0,1 (uncorrected,corrected) bb-purity correction application


        Please verify that: gROOT->GetListOfFiles()->GetSize()>1 otherwise signal template is not added



**/


bool selectSignalTemplates(int _novl=0)
{

std::string tSel("Trig1");

if (!arrayOfHistos)
{arrayOfHistos = new TList();
arrayOfHistos->SetOwner(kTRUE); ///To auto-delete all content of the container
}
if (!arrayOfNames) {
arrayOfNames = new TList();
arrayOfNames->SetOwner(kTRUE); ///To auto-delete all content of the container
}

///Read histograms from file
TFile * file1=0;
 if (gROOT->GetListOfFiles()->GetSize()>1) file1=(TFile*) gROOT->GetListOfFiles()->At(1);


///Parameters of reading
  const int novl = 2;

  TH1F* mjjTemplate[novl];
  const std::string sovl[novl] = { "TCHPT", "TCHP6" };

///Reading histo MC (templates) for signal

_isOk=true; /// assume that all histograms are available

for (int iovl=0; iovl<novl; ++iovl)
{
///mjjbtTCHPTTrig1
if (file1!=0) mjjTemplate[iovl] = (TH1F*) file1->Get(Form("mjjbt%s%s",sovl[iovl].c_str(),tSel.c_str()));
else
{
 mjjTemplate[iovl]  =0 ;

}


}

if (mjjTemplate[_novl] !=0)
        {
                mjjTemplate[_novl]->SetDirectory(gROOT); ///To make histograms independent on _file1;
                arrayOfHistos->Add(mjjTemplate[_novl]);


///Form Name of the template:
         arrayOfNames->Add(new TObjString(Form("signal_%s%s",sovl[_novl].c_str(),tSel.c_str())));
        }
else {
std::cout << "Histogram not found: " << Form("mjjbt%s%s",sovl[_novl].c_str(),tSel.c_str())<< std::endl;
_isOk=false;
}


if (!_isOk ) {std::cout<<"There are problems to read all needed signal templates from the file, error"<<std::endl;
return false;}
_isOk=false; /// reset to ini value


return true;
}



/**
        Parameters of the function:
        _novl = 0,1 (TCHPT, TCHP6), Btag operation point of mass templates
        _corr= 0,1 (uncorrected,corrected) bb-purity correction application


**/


bool selectTemplates(int _novl=0, int _corr=0)
{

std::string tSel("Trig0");
if (!arrayOfHistos) {
arrayOfHistos = new TList();
arrayOfHistos->SetOwner(kTRUE); ///To auto-delete all content of the container
}
if (!arrayOfNames) {
arrayOfNames = new TList();
arrayOfNames->SetOwner(kTRUE); ///To auto-delete all content of the container
}
///Read histograms from file
TFile * file0=0;
 if (gROOT->GetListOfFiles()->GetSize()>0) file0=(TFile*) gROOT->GetListOfFiles()->At(0);

///Parameters of reading
  const int ncateg = 3;
  const int ncorr=2;
  const int ntpat=1;
  const int nfc = 3;
  const int novl = 2;

  TH1F* mjjTemplate[nfc][novl][ncateg][ncorr][ntpat];
//  const std::string tFlav[nfc] = {"Uds", "C", "B"};
  const std::string sFlav[nfc] = {"Q", "C", "B"};

  const std::string sfc[nfc] = { "q", "c", "b" };
  const std::string sovl[novl] = { "TCHPT", "TCHP6" };

///Reading Data

if (file0==0) { cout<<"Can't read data, error"<<endl;return false;}

 TH1F* mjjdata[novl];
  for (int iovl=0; iovl<novl; ++iovl)
  if (file0!=0)  mjjdata[iovl] = (TH1F*) file0->Get(Form("mjjbt%s%s",sovl[iovl].c_str(),tSel.c_str()));
  else mjjdata[iovl] =0;
///Putting Data into arrayOfHistos
for (int iovl=0; iovl<novl; ++iovl)
if (mjjdata[iovl]!=0)
{
        mjjdata[iovl]->SetDirectory(gROOT); ///To make histograms independent on _file0;
        if (iovl==_novl)  {
                arrayOfHistos->Add(mjjdata[iovl]);
                arrayOfNames->Add(new TObjString(Form("data_%s%s",sovl[iovl].c_str(),tSel.c_str())));
                _isOk=true;
        }

}

if (!_isOk) {std::cout<<"The data mass spectrum of "<<sovl[_novl]<<" operation point was not found, error"<<std::endl;
return false;
}
_isOk=false; ///reset to ini value


///Reading histo MC (templates)

_isOk=true; /// assume that all histograms are available
  for (int ifc=0; ifc<nfc; ++ifc)
    for (int iovl=0; iovl<novl; ++iovl)
      for (int icateg=0; icateg<ncateg; ++icateg)
        for (int icorr=0; icorr<ncorr; ++icorr)
          for (int itpat=0; itpat<ntpat; ++itpat)
{

if (file0!=0) 
mjjTemplate[ifc][iovl][icateg][icorr][itpat] = (TH1F*) file0->Get(Form("mjjTemp_%s_%s_Cat%dCorr%dTpat%d%s",sfc[ifc].c_str(),sovl[iovl].c_str(),icateg,icorr,itpat,tSel.c_str()));
else
{
 mjjTemplate[ifc][iovl][icateg][icorr][itpat]  =0 ;
}


}



for (unsigned int i=0;i<sizeof(flavors)/sizeof(int);i++) {
        if (mjjTemplate[flavors[i]][_novl][positions[i]][_corr][0]!=0)
        {
         mjjTemplate[flavors[i]][_novl][positions[i]][_corr][0]->SetDirectory(gROOT); ///To make histograms independent on _file0;
        arrayOfHistos->Add(mjjTemplate[flavors[i]][_novl][positions[i]][_corr][0]);

///Form Name of the template:
        TString _tmpl = "bbb";
        _tmpl.Replace(positions[i],1,sFlav[flavors[i]].c_str());

         arrayOfNames->Add(new TObjString(Form("%s_%s%s",_tmpl.Data(),sovl[_novl].c_str(),tSel.c_str())));

        }

        else {
 std::cout << "Histogram not found: " << Form("mjjTemp_%s_%s_Cat%dCorr%dTpat%d%s",sfc[flavors[i]].c_str(),sovl[_novl].c_str(),positions[i],_corr,0,tSel.c_str())<< std::endl;
_isOk=false;
        }

}
if (!_isOk ) {std::cout<<"There are problems to read all needed templates from the file, error"<<std::endl;
return false;}
_isOk=false; /// reset to ini value


return true;
}




void myFit(int _novl=0, int _corr=0, int _scenario=0,TString _merging="")
{

#ifdef doLogFile
freopen(LOGFILE,"w",stdout);
dup2(fileno(stdout), fileno(stderr));
#endif


TFile  *file=TFile::Open("fitresult.root","UPDATE");

///Read all needed stuff from the _file0
if (!selectTemplates(_novl,_corr) ) {std::cout<<"Something wrong with reading bkg and data, return"<<std::endl; return;}

if (!selectSignalTemplates(_novl) ) {std::cout<<"Something wrong with reading signal"<<std::endl; }

if (arrayOfHistos && !arrayOfHistos->IsEmpty()) dataHist=(TH1*)arrayOfHistos->At(0);

if (dataHist!=0)
{
///To solve the problem zero statistics in bin
///for (int kk=1;kk<=bins;kk++)  if (dataHist->GetBinError(kk)==0 && dataHist->GetBinContent(kk)==0) dataHist->SetBinError(kk,1);
        minVal=dataHist->GetXaxis()->GetXmin();
        maxVal=dataHist->GetXaxis()->GetXmax();
        bins = dataHist->GetNbinsX();
        dataYield=dataHist->Integral();

        ///Validation
        std::cout<<"Validation of mass data spectrum: "<<std::endl;
        std::cout<<"MinVal = "<<minVal<<std::endl;
        std::cout<<"MaxVal = "<<maxVal<<std::endl;
        std::cout<<"Number of bins = "<<bins<<std::endl;
        std::cout<<"Number of all events (data) = "<<dataYield<<std::endl;

/*
double sum=0;
for (int kk=1;kk<=bins;kk++) {
if (dataHist->GetBinError(kk)==0 && dataHist->GetBinContent(kk)==0)dataHist->SetBinError(kk,1);
cout<<"Bin # "<<kk<<" val ="<<dataHist->GetBinContent(kk)<<" error= "<<dataHist->GetBinError(kk)<<endl;
sum+=dataHist->GetBinContent(kk);

}

cout<<"sum = "<<sum<<endl;
*/


///Here possible merging
if (_merging.Length()>0) {
ifstream f(_merging);

///Parsing of file:
TString  str;

///Read first line
str.ReadLine(f);


while (str.Length()>0 )
{
TObjArray* tokens=str.Tokenize(TString(' '));

if (tokens->GetEntries()>1 && !str.Contains("#")) {



TString mother = ((TObjString*) tokens->At(0))->GetString();
TList * daughters = new TList();
daughters->SetOwner(kTRUE); ///To auto-delete all content of the container

for (int dd=1;dd<tokens->GetEntries();dd++)     daughters->Add(new TObjString(tokens->At(dd)->GetName()));
if ( !templateMerge(mother, daughters, 0)) cout<<"Something wrong with merging!!!"<<endl;
//if ( !templateMerge(mother, daughters, 1)) cout<<"Something wrong with merging!!!"<<endl;

daughters->Clear();
delete daughters;
delete tokens;
} /// token
str.ReadLine(f);
}
} ///if _merging




} else {cout<<"dataHist is empty, return"<<endl; return;}



///Common arrays for fitter

const Int_t numPars =20; /// maximal num of params
//if (arrayOfNames && arrayOfNames->GetSize()>1 ) numPars=  arrayOfNames->GetSize() - 1;
static Double_t vstart[numPars]; 
static Double_t step[numPars] ;

static Double_t yields[numPars]; /// container to catch final contribution

Double_t start_val = 0e0;
if (arrayOfNames->GetSize()>1) start_val = dataYield/(arrayOfNames->GetSize()-1);


for (int i=1; i<arrayOfNames->GetSize();i++)
{
TH1 * tempHis = (TH1*)arrayOfHistos->At(i);
TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();

if (!tempHis) continue;


Double_t integral = tempHis->Integral();
if (integral!=0) tempHis->Scale(1e0/integral);
step[i-1]=0.1;
vstart[i-1]=start_val;
//vstart[i-1]=1e0;

}



/////TMinuit 

if ( _scenario == 0 )
{

TMinuit *ptMinuit = new TMinuit(arrayOfNames->GetSize()-1  );  //initialize TMinuit with a maximum of GetSize() - 1


  //  select verbose level:
  //    default :     (58 lines in this test)
  //    -1 : minimum  ( 4 lines in this test)
  //     0 : low      (31 lines)
  //     1 : medium   (61 lines)
  //     2 : high     (89 lines)
  //     3 : maximum (199 lines in this test)

  ptMinuit->SetPrintLevel();


// set the user function that calculates chi_square (the value to minimize)
  ptMinuit->SetFCN(chi_square);

Double_t arglist[10];
  Int_t ierflg = 0;


///Chi2 test error
arglist[0] = 1;
  ptMinuit->mnexcm("SET ERR", arglist ,1,ierflg);


for (int i=1; i<arrayOfNames->GetSize();i++) {

if (i!=2) ptMinuit->mnparm(i-1,Form("yield_%s",arrayOfNames->At(i)->GetName()),vstart[i-1],step[i-1],0.,dataYield,ierflg);
if (i==2) ptMinuit->mnparm(i-1,Form("yield_%s",arrayOfNames->At(i)->GetName()),vstart[i-1],step[i-1],0.,dataYield,ierflg);

}

///Simple minimum search
///Set the strategy 1 (if you want as much reliable as possible then use =2)
 arglist[0] = 2;
  ptMinuit->mnexcm("SET STRategy", arglist ,1,ierflg);




  // Now ready for minimization step

///May improve the global minimum finding
  arglist[0] = 500; ///how many iteration steps at all
  arglist[1] = 2; ///hypercuber in 5 sigma maximal deltaChi2=5^2=25
  ptMinuit->mnexcm("SEEK", arglist ,2,ierflg);



#ifdef doHESSE

/// It might help migrad to work out the error
  arglist[0] = 5000; /// max calls
  ptMinuit->mnexcm("HESse", arglist ,1,ierflg);

///Print correlation matrix

cout<<"Correlation Table 1"<<"\n";
ptMinuit->mnmatu(1);
#endif

double tot=0;


#ifdef doMIGRAD

  arglist[0] = 5000;
  arglist[1] = 1; ///tolerance
  ptMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

  // Print results
  cout << "\n +MYOUTPUT Print results from MIGRAD \n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
  double fParamVal;
  double fParamErr;
  ptMinuit->GetParameter(i-1,fParamVal,fParamErr);
  std::cout <<Form("yield_%s",arrayOfNames->At(i)->GetName()) <<"=" << fParamVal <<"+-"<<fParamErr <<"\n";

yields[i-1]=fParamVal; ///for final plotting 
tot+=fParamVal;
}

cout<<"\n"<<"\n";

for (int i=1; i<arrayOfNames->GetSize();i++)
{
  double fParamVal;
  double fParamErr;
  ptMinuit->GetParameter(i-1,fParamVal,fParamErr);
  std::cout <<Form("ratio yield_%s",arrayOfNames->At(i)->GetName()) <<"=" << fParamVal/tot <<"+-"<<fParamErr/tot <<"\n";
}

cout<<"\n -MYOUTPUT"<<"\n";

#endif

#ifdef doHESSE


/// It might help migrad to work out the error
  arglist[0] = 5000; /// max calls
  ptMinuit->mnexcm("HESse", arglist ,1,ierflg);

///Print correlation matrix

cout<<"Print result from HESSE"<<"\n";
ptMinuit->mnmatu(1);

#endif

#ifdef doMINOS

arglist[0] = 5000; ///max calls =500
ptMinuit->mnexcm("MINOs",arglist,1,ierflg);
tot=0;
  // Print results
  cout << "\n +MYOUTPUT Print results from MINOS\n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
  double fParamVal;
  double fParamErr;
  ptMinuit->GetParameter(i-1,fParamVal,fParamErr);
  std::cout <<Form("yield_%s",arrayOfNames->At(i)->GetName()) <<"=" << fParamVal <<"+-"<<fParamErr <<"\n";
tot+=fParamVal;
}

cout<<"\n"<<"\n";

for (int i=1; i<arrayOfNames->GetSize();i++)
{
  double fParamVal;
    double fParamErr;
  ptMinuit->GetParameter(i-1,fParamVal,fParamErr);
  std::cout <<Form("ratio yield_%s",arrayOfNames->At(i)->GetName()) <<"=" << fParamVal/tot <<"+-"<<fParamErr/tot <<"\n";
}


cout<<"\n -MYOUTPUT"<<"\n";

#endif


// if you want to access to these parameters, use:
  Double_t amin,edm,errdef;
  Int_t nvpar,nparx,icstat;
  ptMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
  //void mnstat(Double_t &fmin, Double_t &fedm, Double_t &errdef, Int_t &npari, Int_t &nparx, Int_t &istat) 
  //*-*-*-*-*Returns concerning the current status of the minimization*-*-*-*-*
  //*-*      =========================================================
  //*-*       User-called
  //*-*          Namely, it returns:
  //*-*        FMIN: the best function value found so far
  //*-*        FEDM: the estimated vertical distance remaining to minimum
  //*-*        ERRDEF: the value of UP defining parameter uncertainties
  //*-*        NPARI: the number of currently variable parameters
  //*-*        NPARX: the highest (external) parameter number defined by user
  //*-*        ISTAT: a status integer indicating how good is the covariance
  //*-*           matrix:  0= not calculated at all
  //*-*                    1= approximation only, not accurate
  //*-*                    2= full matrix, but forced positive-definite
  //*-*                    3= full accurate covariance matrix
  //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
 std::cout << "\n";
  std::cout << " Minimum chi square = " << amin << "\n";
  std::cout << " Estimated vert. distance to min. = " << edm << "\n";
  std::cout << " Number of variable parameters = " << nvpar << "\n";
  std::cout << " Highest number of parameters defined by user = " << nparx << "\n";
  std::cout << " Status of covariance matrix = " << icstat << "\n";

  cout << "\n";
  ptMinuit->mnprin(3,amin);
  //*-*-*-*Prints the values of the parameters at the time of the call*-*-*-*-*
  //*-*    ===========================================================
  //*-*        also prints other relevant information such as function value,
  //*-*        estimated distance to minimum, parameter errors, step sizes.
  //*-*
  //*-*         According to the value of IKODE, the printout is:
  //*-*    IKODE=INKODE= 0    only info about function value
  //*-*                  1    parameter values, errors, limits
  //*-*                  2    values, errors, step sizes, internal values
  //*-*                  3    values, errors, step sizes, first derivs.
  //*-*                  4    values, parabolic errors, MINOS errors
  //*-*    when INKODE=5, MNPRIN chooses IKODE=1,2, or 3, according to ISW(2)
  //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

Int_t npars = ptMinuit->GetNumPars();
Double_t _vals[100];
Double_t _errs[100];
for (int ii=0;ii<npars;ii++)
         ptMinuit->GetParameter(ii, _vals[ii],_errs[ii]);




for (int i=1; i<arrayOfNames->GetSize();i++)
{
  double fParamVal;
  double fParamErr;
//  ptMinuit->GetParameter(i-1,fParamVal,fParamErr);
  std::cout <<Form("yield_%s",arrayOfNames->At(i)->GetName()) <<"=" << _vals[i-1] <<"+-"<<_errs[i-1] <<"\n";

fParamVal =  _vals[i-1];
fParamErr = _errs[i-1];

  arglist[0] = i; ///parameter 1

arglist[1] = NUMSCAN; ///100 points  is maximum!
///Do instead of increasing number of NUMSCANS, narrowing the range of scanning

Int_t numpar=2;

#ifdef NARROWSCAN

//  double fParamVal ;
//  double fParamErr ;
//  ptMinuit->GetParameter(i-1,fParamVal,fParamErr);



#ifdef NARROWRANGE

arglist[2] = fParamVal - NARROWRANGE*fParamErr; /// from bestValue -NARROWRANGE*sigma
arglist[3] = fParamVal + NARROWRANGE*fParamErr; /// to bestValue+NARROWRANGE*sigma
	
	

#else
arglist[2] = fParamVal - 5*fParamErr; /// from bestValue -5sigma
arglist[3] = fParamVal + 5*fParamErr; /// to bestValue+5sigma

#endif

numpar=4;
#endif


///Test Chi2 scan at different point

ptMinuit->mnexcm("SCAN",arglist,numpar,ierflg);

TCanvas *c1 = new TCanvas(Form("CHI2_%s",arrayOfNames->At(i)->GetName()),Form("CHI2_%s",arrayOfNames->At(i)->GetName()));

c1->cd();

///ptMinuit->mnscan(); //this is an internal function
TGraph *gr = (TGraph*)ptMinuit->GetPlot(); ///Get Chi2 plot scan
if (gr){

gr->SetName(Form("gr_chi2_%d_%d_%d_%s",_novl,_corr,_scenario,arrayOfNames->At(i)->GetName()));
gr->SetTitle(Form("chi2_%s",arrayOfNames->At(i)->GetName()));



Double_t el11;
Double_t el21;
Double_t yel11;
Double_t yel21;
Double_t xmin1;
Double_t ymin1;

Double_t el1;
Double_t el2;
Double_t yel1;
Double_t yel2;
Double_t xmin;
Double_t ymin;

cout<<"\n +MYOUTPUT"<<"\n";

TGraph * gr1=0;
gr1 = TMinuitAnalysis(ptMinuit,i-1, fParamVal - 5*fParamErr,  fParamVal + 5*fParamErr,0.,dataYield,el11,el21,yel11,yel21,xmin1,ymin1);

cout<<"Estimated el1 = "<<el11<<"  el2 = "<<el21<<"\n";
cout<<"Estimated yel1 = "<<yel11<<"  yel2 = "<<yel21<<"\n";
cout<<"Estimated xmin = "<<xmin1<<"  ymin = "<<ymin1<<"\n";
cout<<"Estimated xmin -el1 = "<<xmin1-el11<<"\n";
cout<<"Estimated el2 - xmin = "<<el21-xmin1<<"\n";


TGraphAnalysis(gr,el1,el2,yel1,yel2,xmin,ymin);


cout<<"Estimated el1 = "<<el1<<"  el2 = "<<el2<<"\n";
cout<<"Estimated yel1 = "<<yel1<<"  yel2 = "<<yel2<<"\n";
cout<<"Estimated xmin = "<<xmin<<"  ymin = "<<ymin<<"\n";
cout<<"Estimated xmin -el1 = "<<xmin-el1<<"\n";
cout<<"Estimated el2 - xmin = "<<el2-xmin<<"\n";


cout<<"\n -MYOUTPUT"<<"\n";

//gr->SetMarkerColor(kBlue);
//gr->SetMarkerSize(1.5);
//gr->SetMarkerStyle(21);
//gr->SetLineColor(kBlue);

//gr->Draw("a*");
TGraphLinePlot(gr,kBlue, el1,el2,xmin);


//cout<<"gr1 = "<<gr1<<"\n";
if (gr1) {

gr1->SetMarkerColor(kRed);
gr1->SetMarkerSize(1.1);
gr1->SetMarkerStyle(22);
//gr1->SetLineColor(kRed);
//gr1->Draw("P same");
TGraphLinePlot(gr1,kRed, el11,el21,xmin1,gr->GetYaxis(),"P same",0.4);

}


TF1 * parabola=0;
TGraphGetParabola(gr, parabola);

cout<<"parabola="<<parabola<<"\n";

if (parabola)
{
parabola->SetLineColor(kBlue);
parabola->Draw("same");
}

TF1 * parabola1=0;
TGraphGetParabola(gr1, parabola1);
if (parabola1)
{
parabola1->SetLineColor(kRed);
//parabola1->Draw("same");
}

c1->Update();
gPad->Update();


gr->Write();

}


for (int j=i+1; j<arrayOfNames->GetSize();j++)
{
cout<<Form("Under consideration CONT_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName())<<"\n";

TCanvas *c2 = new TCanvas(Form("CONT_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName()),Form("CONT_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName()));
c2->cd();
TGraph *gr2 = (TGraph*)ptMinuit->Contour(60, i-1,j-1); ///Get Chi2 plot scan
TH2D *Hist2D =0;


if (gr2){
gr2->SetName(Form("gr2_CONT_%d_%d_%d_%s_%s",_novl,_corr,_scenario,arrayOfNames->At(i)->GetName(), arrayOfNames->At(j)->GetName()));
gr2->SetTitle(Form("CONT0_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName()));

Double_t a;
Double_t b;
Double_t xmin;
Double_t ymin;
double fParamErr;

xmin=_vals[i-1];
ymin=_vals[j-1];


//  ptMinuit->GetParameter(i-1,xmin,fParamErr);
//  ptMinuit->GetParameter(j-1,ymin,fParamErr);


TEllipseAnalysis(gr2,a,b,xmin,ymin);

gr2->Draw("a*");
gr2->Write();



} else
{

/**
	Minuit was not able to find 4 points lets do it ourselves

**/


//TGraph* gr3 = TMinuitCreateContour( ptMinuit,i-1,j-1,arrayOfNames,2,1,0.,dataYield );
//if (gr3) {
//TGraphEllipsePlot(gr3,kBlue,kRed,_vals[i-1],_vals[j-1]);
//} //if


Double_t  _vals0[100];
Double_t  _errs0[100];

for (int ii=0;ii<npars;ii++){
	_vals0[ii]=_vals[ii];
}

Double_t _fval;

ptMinuit->Eval(npars, 0, _fval, _vals,ierflg);


const Int_t nRun=10000;
const Int_t nSigma=2;
const Int_t nBin2D=30;
Double_t minX=0;
Double_t maxX=0;
Double_t minY=0;
Double_t maxY=0;
Double_t contours[nSigma-1]; /// 1,2,3 ...  nSigma-1 sigma contours

for (int k=0;k<nSigma-1;k++)
contours[k] = _fval+(k+1)*(k+1);

minX=_vals[i-1]-nSigma*_errs[i-1];
maxX=_vals[i-1]+nSigma*_errs[i-1];


if (minX<0.) minX=0.;
if (maxX>dataYield) maxX=dataYield;

minY=_vals[j-1]-nSigma*_errs[j-1];
maxY=_vals[j-1]+nSigma*_errs[j-1];


if (minY<0.) minY=0.;
if (maxY>dataYield) maxY=dataYield;




Hist2D = new TH2D(Form("H2CONT_%s_%s",arrayOfNames->At(i)->GetName(),
arrayOfNames->At(j)->GetName()),Form("H2CONT_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName()),nBin2D,minX,maxX,nBin2D,minY,maxY);
Hist2D->SetDirectory(gROOT);

//Hist2D->SetContour(nSigma-1, contours);
//Hist2D->SetContour(1, contours);

///New 
arglist[0] =  i;
ptMinuit->mnexcm("FIX", arglist ,1,ierflg);
arglist[0] =  j;
ptMinuit->mnexcm("FIX", arglist ,1,ierflg);


for (int k=0;k<nRun;k++)
{

Double_t _X = gRandom->Uniform(minX,maxX);
Double_t _Y = gRandom->Uniform(minY,maxY);

///New ! Perform finding optimal 


//_vals0[i-1]=_X;
//_vals0[j-1]=_Y;

 arglist[0] = i;
 arglist[1] =  _X;
 ptMinuit->mnexcm("SET PAR", arglist ,2,ierflg);
 arglist[0] = j;
 arglist[1] =  _Y;
 ptMinuit->mnexcm("SET PAR", arglist ,2,ierflg);


 arglist[0] = 500;
  arglist[1] = 1; ///tolerance
  ptMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

for (int kk=0;kk<npars;kk++)
        ptMinuit->GetParameter(kk, _vals0[kk],_errs0[kk]);

ptMinuit->Eval(npars, 0, _fval, _vals0,ierflg);

for (int iii=0;iii<npars;iii++)
ptMinuit->DefineParameter(iii,Form("yield_%s",arrayOfNames->At(iii+1)->GetName()),_vals[iii],_errs[iii],0.,dataYield);








//Hist2D->Fill(_X,_Y,_fval);

if (_fval<=contours[0]) {
 Hist2D->Fill(_X,_Y);
}

} /// for nRun

///Restore all initial vals 
arglist[0] =  i;
ptMinuit->mnexcm("RELease", arglist ,1,ierflg);
arglist[0] =  j;
ptMinuit->mnexcm("RELease", arglist ,1,ierflg);


//Hist2D->Draw("cont2,list");
Hist2D->Draw("SCAT");
c2->Update();
gPad->Update();

///Calculate TGraph for analysis and plottings

Double_t DownLimitX[nBin2D];
Double_t DownLimitY[nBin2D];
Double_t UpLimitX[nBin2D];
Double_t UpLimitY[nBin2D];
Int_t numPointsDown=0;
Int_t numPointsUp=0;
Int_t row=0;

for (int k1=1;k1<=nBin2D;k1++)
{
Int_t numInternalPoints=0;
for (int k2=1;k2<=nBin2D;k2++)
{

Double_t _cont = Hist2D->GetBinContent(k1,k2);
Double_t xcnt;
Double_t ycnt;

if (_cont!=0) {
numInternalPoints++;
TAxis *xaxis = Hist2D->GetXaxis();
TAxis *yaxis = Hist2D->GetYaxis();
xcnt = xaxis->GetBinCenter(k1);
ycnt = yaxis->GetBinCenter(k2);


if (numInternalPoints==1) 
{
	DownLimitX[numPointsDown]=xcnt;
	DownLimitY[numPointsDown]=ycnt;
	numPointsDown++;
} 

if (numInternalPoints>1)
{
        UpLimitX[numPointsUp]=xcnt;
        UpLimitY[numPointsUp]=ycnt;

}


}


if (k2==nBin2D && numInternalPoints>1) numPointsUp++;

} /// for k2
} /// for k1


TGraph *gr3 =0;
if (numPointsDown+numPointsUp>0) {
gr3= new TGraph(numPointsDown+numPointsUp);
for (int k1=0; k1<numPointsDown;k1++)
gr3->SetPoint(k1,DownLimitX[k1],DownLimitY[k1]);
for (int k2=0;k2<numPointsUp;k2++)
gr3->SetPoint(numPointsDown+k2,UpLimitX[k2],UpLimitY[k2]);

gr3->SetName(Form("gr3_CONT_%d_%d_%d_%s_%s",_novl,_corr,_scenario,arrayOfNames->At(i)->GetName(), arrayOfNames->At(j)->GetName()));
gr3->SetTitle(Form("CONT_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName()));

TCanvas *c3 = new TCanvas(Form("CONT2_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName()),Form("CONT2_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName()));


Double_t a;
Double_t b;
Double_t xmin;
Double_t ymin;

xmin=_vals[i-1];
ymin=_vals[j-1];

TGraphEllipsePlot(gr3,kBlue,kRed,xmin,ymin);


//gr3->Draw("a*");
gr3->Write();


//TEllipseAnalysis(gr3,a,b,xmin,ymin);



c3->Update();
gPad->Update();
}



/*
// Get Contours
   TObjArray *conts = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
cout<<"conts="<<conts<<"\n";
TCanvas *c3 = new TCanvas(Form("CONT2_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName()),Form("CONT2_%s_%s",arrayOfNames->At(i)->GetName(),arrayOfNames->At(j)->GetName()));
c3->cd();


for (int ic=0; ic<conts->GetEntries() && ic<1; ic++) {    ///plot 1st grap
    TList* graphlist = (TList*)conts->At(ic);   
   cout<<"Number of graphs = "<<graphlist->GetEntries()<<endl;
    for (int ig=0; ig<graphlist->GetEntries(); ig++) {     

      TGraph* gr3 = (TGraph*)graphlist->At(ig);
	cout<<"Number of points = "<<gr->GetN()<<endl;
	if (ig>0)      gr3->Draw("a* same");
	else   gr3->Draw("a*");

    }
  }

c3->Update();
gPad->Update();



///To restore

///Allow printings
ptMinuit->mnexcm("SET WARnings", arglist ,0,ierflg);
ptMinuit->SetPrintLevel();

for (int ii=0;ii<npars;ii++)
{
// arglist[0] =  ii+1;
// arglist[1] =  _vals[ii];
// ptMinuit->mnexcm("SET PAR", arglist ,2,ierflg);

ptMinuit->DefineParameter(ii,Form("yield_%s",arrayOfNames->At(ii+1)->GetName()),_vals[ii],_errs[ii],0.,dataYield);


}

 /// Restore HESSE
  arglist[0] = 3000;
  arglist[1] = 1; ///tolerance
  ptMinuit->mnexcm("MINOS", arglist ,2,ierflg);
*/


} ///else

} ///for j
} ///for i



///Extract Correlation matrix 

#ifdef doHESSE


/// It might help migrad to work out the error
  arglist[0] = 5000; /// max calls
  ptMinuit->mnexcm("HESse", arglist ,1,ierflg);

///Print correlation matrix

cout<<"Print result from HESSE"<<"\n";
ptMinuit->mnmatu(1);

#endif



//Double_t _corrMaxtrix[arrayOfNames->GetSize()-1][arrayOfNames->GetSize()-1];
//Double_t _corrMaxtrix[9][9];

Int_t npars = ptMinuit->GetNumPars();
//Double_t *covar = new Double_t[npars*npars]; 

TMatrixDSym _corrMaxtrix(npars); 


cout<<"\n +MYOUTPUT"<<"\n";
cout<<"External correlation Matrix from TMinuit"<<"\n";

//ptMinuit->mnemat(&_corrMaxtrix[0][0],arrayOfNames->GetSize()-1 );
//ptMinuit->mnemat(covar,npars );
ptMinuit->mnemat(_corrMaxtrix.GetMatrixArray(),npars );


//for (int i=1;i<arrayOfNames->GetSize();i++)
//cout<<"a"<<i<<"0 = "<<covar[i][0]<<" a"<<i<<"1 = "<<covar[i][1]<<" a"<<i<<"2 = "<<covar[i][2]<<"\n";
//cout<<"a"<<i<<"0 = "<<_corrMaxtrix[i-i][0]<<" a"<<i<<"1 = "<<_corrMaxtrix[i-i][1]<<" a"<<i<<"2 = "<<_corrMaxtrix[i-i][2]<<"\n";
for (int i=0;i<npars;i++){
for (int j=0;j<npars;j++)
cout<<"a"<<i<<j<<"="<<_corrMaxtrix(i,j)<<"   ";

cout<<"\n";
}
///Internal error matrix
/*
TVirtualFitter * fitter =  TVirtualFitter::GetFitter();
for (int i=1;i<=npars;i++)
for (int j=1;j<=npars;j++)
{

 double corr_ij = fitter->GetCovarianceMatrixElement(i,j) / sqrt( fitter->GetCovarianceMatrixElement(i,i) *
                                                   fitter->GetCovarianceMatrixElement(j,j) ) ;

cout<< "cor["<<i<<","<<j<<"]"<<corr_ij<<"\n";
}
*/

///Fix parameter 1

ptMinuit->FixParameter(0);
///Do HESSE again
/// It might help migrad to work out the error
  arglist[0] = 5000; /// max calls
  ptMinuit->mnexcm("HESse", arglist ,1,ierflg);

///Print correlation matrix

cout<<"Print result from HESSE"<<"\n";
ptMinuit->mnmatu(1);



///Do instead of increasing number of NUMSCANS, narrowing the range of scanning



cout<<"\n -MYOUTPUT"<<"\n";



///delete ptMinuit;
} /// if 


////TFitter case
if (_scenario == 1)
{

TFitter* minimizer = new TFitter(arrayOfNames->GetSize()-1  );

// MAKE IT MEDIUM LEVEL VERBOSITY!!

double p1 = 1;
minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);

minimizer->SetFCN(chi_square);

for (int i=1; i<arrayOfNames->GetSize();i++) 
minimizer->SetParameter(i-1,Form("yield_%s",arrayOfNames->At(i)->GetName()),vstart[i-1],1,0.,dataYield);

// Run the migrad minimizer (an extended Powell's method) to improve the
// fit.
minimizer->ExecuteCommand("MIGRAD",0,0);


//Needed to get the function value at the best fit.
double minimum ;
Int_t  npar = arrayOfNames->GetSize() -1; /// all except data
Double_t pars[numPars];
Double_t globalpars[numPars];
Int_t ierflg = 0;


double _totyield=0;

// Print results
  cout << "\n +MYOUTPUT Print results from TFitter\n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
  double fParamVal = minimizer->GetParameter(i-1);
  double fParamErr = minimizer->GetParError(i-1);
  std::cout <<Form("yield_%s",arrayOfNames->At(i)->GetName()) <<"=" << fParamVal <<"+-"<<fParamErr <<"\n";

_totyield+=fParamVal;

pars[i-1] = fParamVal;
globalpars[i-1] = fParamVal;
yields[i-1]=fParamVal; ///for final plotting

}

cout<<"\n -MYOUTPUT"<<"\n";

// Print results
  cout << "\n +MYOUTPUT Print ratios from TFitter\n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
  double fParamVal = minimizer->GetParameter(i-1);
  double fParamErr = minimizer->GetParError(i-1);
  std::cout <<Form("ratio_yield_%s",arrayOfNames->At(i)->GetName()) <<"=" << fParamVal/_totyield <<"+-"<<fParamErr/_totyield <<"\n";
}


cout<<"\n -MYOUTPUT"<<"\n";


///Derived method from MINOS (Chi2 scan)


chi_square(npar,0,minimum,pars,ierflg);

cout<<"\n +MYOUTPUT"<<"\n";


std::cout <<"Global minimum = "<<minimum<<"\n";

cout<<"\n -MYOUTPUT"<<"\n";


const Int_t _numsteps=500;



///No verbosity at all
double p1 = -1;
minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);



///loop over parameters 
for (int i=1; i<arrayOfNames->GetSize();i++)
{

double _xmin=0.;
double _xmax=dataYield;
double _step = (_xmax - _xmin)/_numsteps;
_xmax+=_step;

///needed for estimation of the error
Double_t _valMinus=0;
Double_t _valPlus=0;

Double_t _deltachi2_y[_numsteps+1];
Double_t _deltachi2_x[_numsteps+1];


Int_t _cur_step=0;
// Scan the Chi2 over  parameter to find it's uncertainty
double _parVal;
for (_parVal=_xmin; _parVal<_xmax; _parVal = _parVal + _step) {

minimizer->SetParameter(i-1,Form("yield_%s",arrayOfNames->At(i)->GetName()),_parVal,0,0.,0.);

for (int j=1; j<arrayOfNames->GetSize();j++)
if (i !=j ) minimizer->SetParameter(j-1,Form("yield_%s",arrayOfNames->At(j)->GetName()),vstart[j-1],1,0.,dataYield);

// Run the migrad minimizer (an extended Powell's method) to improve the
// fit.
minimizer->ExecuteCommand("MIGRAD",0,0);

for (int k=1; k<arrayOfNames->GetSize();k++)
{
  double fParamVal = minimizer->GetParameter(k-1);
 pars[k-1] = fParamVal;
}

double t;
chi_square(npar,0,t,pars,ierflg);

_deltachi2_y[_cur_step] = t-minimum;
_deltachi2_x[_cur_step] = _parVal;

///if (abs(t-minimum) <= 1.0) break;

/*
cout<<"Var = "
<< Form("yield_%s",arrayOfNames->At(i)->GetName())
<<" val = "
<<_parVal <<" current step = "<<_cur_step<<" deltaChi2 = "<< t-minimum <<"\n";
*/

_cur_step++;


} /// scan

///Look for vals when (t-minimum) ~ 1.0)

int flag=-1; ///determines 3 different cases: val =0 , when minimum close to the left limit, val=1, when the minimum in the middle of the range
/// val=2, when the minimum at right limit

if (_deltachi2_y[1]>1e0 && _deltachi2_y[_numsteps]>1e0  ) flag=1;
if (_deltachi2_y[1]<1e0 && _deltachi2_y[_numsteps]>1e0  ) flag=0;
if (_deltachi2_y[1]>1e0 && _deltachi2_y[_numsteps]<1e0  ) flag=2;

_valMinus=-1e30;
_valPlus=-1e30;

//cout<<"Case "<<flag<<endl;

for (int m=1;m<_numsteps;m++) {


if (_deltachi2_y[m] < 1e0 && flag==0 ) {
_valMinus = _deltachi2_x[m];
///cout<<"I found new min "  << " at - "<<_valMinus<<"\n";
//break;
}

if (_deltachi2_y[m] < 1e0 && flag==2 && _valPlus<0 ) {
_valPlus = _deltachi2_x[m];
///cout<<"I found new min "  << " at + "<<_valPlus<<"\n";
//break;
}

if (flag==1) {
if (_deltachi2_y[m] < 1e0 && _valMinus<0) _valMinus=_deltachi2_x[m];
if (_deltachi2_y[m] < 1e0 ) _valPlus=_deltachi2_x[m];
}

}

/*
double min=1e30;
for (int m=1;m<_numsteps;m++)
if (abs(_deltachi2_y[m] - 1e0)< min ) {
min=abs(_deltachi2_y[m] - 1e0);
_valPlus = _deltachi2_x[m];
cout<<"I found new min at + "<<_valPlus<<"\n";

}
*/

double tmp=_valPlus;
if (_valMinus<0 && _valPlus>0) _valMinus=_valPlus;
if (_valMinus>0 && _valPlus<0) _valPlus=_valMinus;

if (_valMinus>_valPlus) {_valPlus=_valMinus; _valMinus=tmp;};
  std::cout <<"\n\n +MYOUTPUT "<<Form("yield_%s",arrayOfNames->At(i)->GetName()) <<"= "<<globalpars[i-1] <<" + "<< abs(_valPlus - globalpars[i-1])<<" - "<<abs(globalpars[i-1] - _valMinus) <<"\n\n -MYOUTPUT";


} /// loop over parameters

///cout<<"Stop"<<endl;
delete minimizer;
}



////TFitter case
if (_scenario == 2)
{

TFitter* minimizer = new TFitter(arrayOfNames->GetSize()-1  );

// MAKE IT MEDIUM LEVEL VERBOSITY!!

double p1 = 1;
minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);

minimizer->SetFCN(chi_square);

for (int i=1; i<arrayOfNames->GetSize();i++) 
minimizer->SetParameter(i-1,Form("yield_%s",arrayOfNames->At(i)->GetName()),vstart[i-1],1,0.,dataYield);

// Run the migrad minimizer (an extended Powell's method) to improve the
// fit.
minimizer->ExecuteCommand("MIGRAD",0,0);


//Needed to get the function value at the best fit.
double minimum ;
Int_t  npar = arrayOfNames->GetSize() -1; /// all except data
Double_t pars[numPars];
Double_t globalpars[numPars];
Int_t ierflg = 0;


double _totyield=0;

// Print results
  cout << "\nPrint results from TFitter\n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
  double fParamVal = minimizer->GetParameter(i-1);
  double fParamErr = minimizer->GetParError(i-1);
  std::cout <<Form("yield_%s",arrayOfNames->At(i)->GetName()) <<"=" << fParamVal <<"+-"<<fParamErr <<"\n";

_totyield+=fParamVal;

pars[i-1] = fParamVal;
globalpars[i-1] = fParamVal;
yields[i-1]=fParamVal; ///for final plotting
}

// Print results
  cout << "\nPrint ratios from TFitter\n";
for (int i=1; i<arrayOfNames->GetSize();i++)
{
  double fParamVal = minimizer->GetParameter(i-1);
  double fParamErr = minimizer->GetParError(i-1);
  std::cout <<Form("ratio_yield_%s",arrayOfNames->At(i)->GetName()) <<"=" << fParamVal/_totyield <<"+-"<<fParamErr/_totyield <<"\n";
}


///Derived method from MINOS (Chi2 scan)


chi_square(npar,0,minimum,pars,ierflg);
std::cout <<"Global minimum = "<<minimum<<"\n";


const Int_t _numsteps=500;



///No verbosity at all
double p1 = -11;
minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);



///loop over parameters 
for (int i=1; i<arrayOfNames->GetSize();i++)
{

double _xmin=0.;
double _xmax=dataYield;
double _step = (_xmax - _xmin)/_numsteps;
_xmax+=_step;

///needed for estimation of the error
Double_t _valMinus=0;
Double_t _valPlus=0;

Double_t _deltachi2_y[_numsteps+1];
Double_t _deltachi2_x[_numsteps+1];


Int_t _cur_step=0;
// Scan the Chi2 over  parameter to find it's uncertainty
double _parVal;
for (_parVal=_xmin; _parVal<_xmax; _parVal = _parVal + _step) {

minimizer->SetParameter(i-1,Form("yield_%s",arrayOfNames->At(i)->GetName()),_parVal,0,0.,0.);

for (int j=1; j<arrayOfNames->GetSize();j++)
///if (i !=j ) minimizer->SetParameter(j-1,Form("yield_%s",arrayOfNames->At(j)->GetName()),vstart[j-1],1,0.,dataYield);
if (i !=j ) minimizer->SetParameter(j-1,Form("yield_%s",arrayOfNames->At(j)->GetName()),globalpars[j-1] ,0,0.,0.);

// Run the migrad minimizer (an extended Powell's method) to improve the
// fit.

minimizer->ExecuteCommand("MIGRAD",0,0);

for (int k=1; k<arrayOfNames->GetSize();k++)
{
  double fParamVal = minimizer->GetParameter(k-1);
 pars[k-1] = fParamVal;
}

double t;
chi_square(npar,0,t,pars,ierflg);

_deltachi2_y[_cur_step] = t-minimum;
_deltachi2_x[_cur_step] = _parVal;

///if (abs(t-minimum) <= 1.0) break;

/*
cout<<"Var = "
<< Form("yield_%s",arrayOfNames->At(i)->GetName())
<<" val = "
<<_parVal <<" current step = "<<_cur_step<<" deltaChi2 = "<< t-minimum <<"\n";
*/

_cur_step++;


} /// scan

///Look for vals when (t-minimum) ~ 1.0)

int flag=-1; ///determines 3 different cases: val =0 , when minimum close to the left limit, val=1, when the minimum in the middle of the range
/// val=2, when the minimum at right limit

if (_deltachi2_y[1]>1e0 && _deltachi2_y[_numsteps]>1e0  ) flag=1;
if (_deltachi2_y[1]<1e0 && _deltachi2_y[_numsteps]>1e0  ) flag=0;
if (_deltachi2_y[1]>1e0 && _deltachi2_y[_numsteps]<1e0  ) flag=2;

_valMinus=-1e30;
_valPlus=-1e30;

//cout<<"Case "<<flag<<endl;

for (int m=1;m<_numsteps;m++) {


if (_deltachi2_y[m] < 1e0 && flag==0 ) {
_valMinus = _deltachi2_x[m];
///cout<<"I found new min "  << " at - "<<_valMinus<<"\n";
//break;
}

if (_deltachi2_y[m] < 1e0 && flag==2 && _valPlus<0 ) {
_valPlus = _deltachi2_x[m];
///cout<<"I found new min "  << " at + "<<_valPlus<<"\n";
//break;
}

if (flag==1) {
if (_deltachi2_y[m] < 1e0 && _valMinus<0) _valMinus=_deltachi2_x[m];
if (_deltachi2_y[m] < 1e0 ) _valPlus=_deltachi2_x[m];
}

}

/*
double min=1e30;
for (int m=1;m<_numsteps;m++)
if (abs(_deltachi2_y[m] - 1e0)< min ) {
min=abs(_deltachi2_y[m] - 1e0);
_valPlus = _deltachi2_x[m];
cout<<"I found new min at + "<<_valPlus<<"\n";

}
*/

double tmp=_valPlus;
if (_valMinus<0 && _valPlus>0) _valMinus=_valPlus;
if (_valMinus>0 && _valPlus<0) _valPlus=_valMinus;

if (_valMinus>_valPlus) {_valPlus=_valMinus; _valMinus=tmp;};
  std::cout <<"\n\n +MYOURPUT "<<Form("yield_%s",arrayOfNames->At(i)->GetName()) <<"= "<<globalpars[i-1] <<" + "<< abs(_valPlus - globalpars[i-1])<<" - "<<abs(globalpars[i-1] - _valMinus) <<"\n\n -MYOUTPUT";



} /// loop over parameters

////cout<<"Stop"<<endl;
///delete minimizer;
} //scenario==2


///Plotting MC stack

///Perform Sorting 
TGraph   _combined;

double _totalMC=0;

for (int i=1; i<arrayOfNames->GetSize();i++)
{
TH1 * tempHis = (TH1*)arrayOfHistos->At(i);
TString _name = ( (TObjString *) arrayOfNames->At(i))->GetString();

if (!tempHis) continue;

_totalMC+=yields[i-1];
tempHis->Scale(yields[i-1]);

_combined.SetPoint(i,i,tempHis->GetMaximum());

}


///Make Sorting
_combined.Sort(&TGraph::CompareY,kFALSE); ///descending on Y
Double_t * templPos = _combined.GetX();


TCanvas * c5 = new TCanvas ("final_plot","final plot");
TLegend * leg = new TLegend(0.4,0.6,0.89,0.89);
THStack *hs = new THStack("hs","final plot");
if (dataHist) leg->AddEntry(dataHist,"data","l");

for (int i=0; i<_combined.GetN();i++)
{
if (templPos[i] == 0 ) continue; ///artifact --> TGraph always adds (0,0) point if it doesn't exist
TString _name = ( (TObjString *) arrayOfNames->At(templPos[i]))->GetString();
TH1 * tempHis = (TH1*)arrayOfHistos->At(templPos[i]);
if (!tempHis) continue;
tempHis->SetLineColor(colors[int(templPos[i])]);
tempHis->SetFillColor(colors[int(templPos[i])]);
tempHis->SetFillStyle(1001);

leg->AddEntry(tempHis,_name,"l");
hs->Add(tempHis);

}

c5->cd();
hs->Draw("hist");
if (dataHist) dataHist->Draw("e0 same");
leg->Draw();

///Save plots
TSeqCollection * canvs = 0;
canvs = gROOT->GetListOfCanvases();


if (canvs)
for (int i=0; i<canvs->GetEntries(); i++)
((TCanvas*)canvs->At(i))->SaveAs(TString(canvs->At(i)->GetName())+TString(".eps"));

file->Close();

cout<<"+MYOUTPUT "<<"\n";
if (dataHist) cout<<"Data : "<<dataHist->Integral()<<endl;
cout<<"Fit : "<<_totalMC<<endl;
cout<<"-MYOUTPUT "<<"\n";


return;

}
